%% ========================================================================
%  MATLAB 演示：从 H(s) 到差分方程的仿真之旅
%  ========================================================================
%  本脚本旨在用最简洁的方式，演示以下核心流程：
%  1. 定义一个连续传递函数 H(s)
%  2. 将其通过Z变换，转换为离散传递函数 H(z) / 差分方程
%  3. 手动实现差分方程，并用它来对一个信号进行滤波
%  4. 将手动实现的结果与MATLAB内置函数进行对比，验证其正确性
% ========================================================================
%【闲鱼】https://m.tb.cn/h.hnuhKvM?tk=sVEc45fMvZk HU926 「我在闲鱼发布了【2025电赛G题题目发挥部分的代码，模拟效果演示，pytho】」点击链接直接打开
%咸鱼ID:代代代鱼
clear; clc; close all;

% --- 为了让中文标题能够正常显示 ---
set(groot, 'defaultAxesFontName', 'Microsoft YaHei');
set(groot, 'defaultTextFontName', 'Microsoft YaHei');


%% 步骤 1: 定义一个我们已知的“物理世界”模型 H(s)
% -------------------------------------------------------------------------
% 假设我们有一个简单的二阶低通滤波器（巴特沃斯滤波器）
% 它的自然谐振频率 wn = 1000 rad/s (约 159 Hz)
% 阻尼比 zeta = 0.707 (临界阻尼)

wn = 1000;
zeta = 0.707;

% 连续传递函数 H(s) = wn^2 / (s^2 + 2*zeta*wn*s + wn^2)
num_s = [wn^2];
den_s = [1, 2*zeta*wn, wn^2];

G_s = tf(num_s, den_s);

disp('--- 1. 连续系统 H(s) ---');
disp('我们从一个已知的连续传递函数开始:');
G_s


%% 步骤 2: 将 H(s) 转换为离散的差分方程
% -------------------------------------------------------------------------
% 我们需要选择一个采样率 Fs。根据奈奎斯特定理，它必须高于系统
% 主要频率的两倍。我们选择一个远高于159Hz的采样率，例如 2000 Hz。
Fs = 2000; % 采样率 (Hz)
Ts = 1/Fs; % 采样周期 (秒)

% 使用双线性变换 (Tustin's method) 将 H(s) 转换为 H(z)
% c2d 函数 (Continuous to Discrete) 是MATLAB的专业工具
G_z = c2d(G_s, Ts, 'tustin');

% 从 G_z 中提取差分方程的 b 和 a 系数
% 在MATLAB中，分子是 B(z)，分母是 A(z)
b = G_z.Numerator{1};
a = G_z.Denominator{1};

disp('--- 2. 离散化 -> 差分方程 ---');
disp(['选定采样率 Fs = ', num2str(Fs), ' Hz']);
disp('使用双线性变换得到 H(z):');
G_z
disp('这直接对应了差分方程的系数：');
fprintf('  分子系数 b = [%.4f, %.4f, %.4f]\n', b(1), b(2), b(3));
fprintf('  分母系数 a = [%.1f, %.4f, %.4f]\n', a(1), a(2), a(3));

% 我们可以写出对应的差分方程 (假设 a(1) = 1)
fprintf('\n对应的差分方程为:\n');
fprintf('  y[n] = b(1)*x[n] + b(2)*x[n-1] + b(3)*x[n-2] - a(2)*y[n-1] - a(3)*y[n-2]\n');


%% 步骤 3: 手动实现差分方程并进行仿真
% -------------------------------------------------------------------------
disp('--- 3. 正在手动执行差分方程... ---');

% 创建一个“不完美”的输入信号，包含基波和高频噪声
t = 0:Ts:1; % 1秒的时间
f_base = 10;  % 10Hz 的基波方波
f_noise = 200; % 200Hz 的高频噪声

u_signal = square(2*pi*f_base*t) + 0.3*sin(2*pi*f_noise*t);

% 初始化输出向量和历史状态变量
y_manual = zeros(size(u_signal));
x_n_1 = 0; % x[n-1]
x_n_2 = 0; % x[n-2]
y_n_1 = 0; % y[n-1]
y_n_2 = 0; % y[n-2]

% --- 这就是单片机中的主循环 (For Loop) ---
for n = 1:length(u_signal)
    % 获取当前输入
    x_n = u_signal(n);
    
    % 应用差分方程 (我们已经归一化 a(1)=1)
    y_n = b(1)*x_n + b(2)*x_n_1 + b(3)*x_n_2 - a(2)*y_n_1 - a(3)*y_n_2;
    
    % 保存当前输出
    y_manual(n) = y_n;
    
    % 更新历史状态，为下一次循环做准备
    % x[n] -> x[n-1], x[n-1] -> x[n-2]
    x_n_2 = x_n_1;
    x_n_1 = x_n;
    
    % y[n] -> y[n-1], y[n-1] -> y[n-2]
    y_n_2 = y_n_1;
    y_n_1 = y_n;
end

disp('手动执行完毕。');

%% 步骤 4: 与MATLAB内置函数对比，验证结果
% -------------------------------------------------------------------------
disp('--- 4. 正在与内置函数 filter() 对比验证... ---');

% MATLAB 的 filter() 函数就是 lfilter 的孪生兄弟，专门用于执行差分方程
y_builtin = filter(b, a, u_signal);

% 计算两种方法之间的误差
error = max(abs(y_manual - y_builtin));
fprintf('手动实现与内置函数的最大误差: %.4e\n', error);
if error < 1e-12
    disp('结果完美匹配！验证成功！');
else
    disp('结果存在差异，请检查代码。');
end

% --- 最终绘图 ---
figure('Name', 'MATLAB 差分方程演示');
plot(t, u_signal, 'b--', 'DisplayName', '输入信号 (带噪声方波)', 'LineWidth', 1);
hold on;
plot(t, y_manual, 'r-', 'DisplayName', '手动实现的差分方程输出', 'LineWidth', 3);
plot(t, y_builtin, 'g:', 'DisplayName', 'MATLAB内置函数输出', 'LineWidth', 2);
hold off;

title('差分方程的实现与验证', 'FontSize', 16);
xlabel('时间 (秒)');
ylabel('幅值');
grid on;
legend('show');
ylim([-1.5, 1.5]);